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I  DISTRIBUTION  STATEMENT  A  I 


Abstract 


This  paper  deals  with  the  stable  c-server  queue  with  renewal 

input.  The  service  time  distributions  may  be  different  for  the 

various  servers.  They  are  however  all  probability  distributions 

of  phase  type.  It  is  shown  that  the  stationary  distribution  of 

the  queue  length  at  arrivals  has  an  exact  geometric  tail  of  rate 

n»  0  <  n  <  1.  It  is  further  shown  that  the  stationary  waiting 

time  distribution  at  arrivals  has  an  exact  exponential  tail  of 

decay  parameter  £  >  0.  The  quantities  r\  and  £  may  be  evaluated 

together  by  an  elementary  algorithm.  For  both  distributions,  the 

multiplicative  constants  which  arise  in  the  asymptotic  forms  may  be 

»  * 

fully  characterized.  These  constants  are  however  difficult  to 
compute  in  general. 
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1.  Introduction 


Very  few  algorithmically  tractable  results  are  known  for 
multi-server  queues,  except  In  the  restrictive  ease  where  the 
service  time  distributions  are  exponential.  In  this  paper,  we 
obtain  results  on  the  tall  behavior  of  the  stationary  distributions 
of  the  queue  length  and  the  waiting  time  for  a  c -server  queue  with 
renewal  input.  The  various  servers  are  allowed  to  be  heterogeneous , 
l.e.  the  service  time  distributions  may  be  different  for  different 
servers.  The  c  service  time  distributions  are  however  required 
to  be  of  phase  type. 

If  Pm  and  W(x)  denote  respectively  the  stationary  proba¬ 
bilities  that  a  customer  arriving  to  the  queue  finds  at  least  m 
customers  in  the  system  and  that  he  has  to  wait  for  a  time  at  most 
x,  then  we  establish  the  asymptotic  formulas 

Pffl  ■  Kqm  +  o(qm),  as  m  -►  «®, 

and 

1  -  W(x)  ■  e“^x  +  o(e"^x),  as  x  •*>  •, 

where  K1  »  n°K,  The  constants  q  and  C,  which  satisfy  0  <  n  <  1 
and  5  >0,  will  be  shown  to  satisfy  a  system  of  equations,  which 
may  easily  be  sqlved  by  elementary  numerical  procedures.  The 
positive  constants  K  and  will  be  fully  characterized,  but 

their  direct  numerical  evaluation  is  seen  to  be  difficult,  in 
general . 
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In  the  course  of  the  discussion,  a  number  of  results  of 
theoretical  Interest  on  the  GI/PH/c  queue  with  heterogeneous 
servers  will  be  obtained.  Because  of  the  huge  dimensions  of  the 
matrices,  which  arise  In  these  results,  their  feasibility  for 
algorithmic  Implementation  Is  limited  to  particular  cases. 

The  proofs  of  the  various  results  in  this  paper  rely  heavily 
on  the  elementary  properties  of  phase  type  distributions  [4,  5,  8] 
and  on  the  theory  of  block-partitioned  stochastic  matrices  with  a 
matrix-geometric  Invariant  vector  [ 6,  7,  8].  Specific  references 
will  be  given  for  each  prior  result,  which  is  used,  but  no  easily 
available  proofs  will  be  repeated.  The  particular  results  for  the 
GI/PH/1  queue  were  established  in  Neuts  C8,  9].  Corresponding 
results  for  the  PH/PH/c  queue  with  identical  servers  are  discussed 
in  Takahashi  [11],  while  iterative  numerical  procedures  for  that 
model  are  proposed  in  Takahashi  and  Takaml  [10]. 

The  formal  description  of  the  model  is  as  follows.  Customers 
arrive  to  a  c-server  system  according  to  a  renewal  process  with 
interarrival  time  distribution  ?(•)  of  finite  mean  X*  .  The 
distribution  P(’)  satisfies  F(0+)  ■  0.  The  service  time  distri¬ 
bution  of  a  customer  may  depend  on  the  server  to  which  he  gains 
access.  Services  by  the  J-th  server,  1  s  J  sc,  have  a  common 
distribution  of  phase  type,  with  irreducible  representation 
[JJ(J),  SCJ)]  where  £(J )  is  a  probability  row-vector  of  dimension 
v(J)  and  S(J)  is  a  square,  stable  matrix  of  order  v(J).  The 
corresponding  vector  S°(J)  is  defined  by  S° (J )  »  -  S(J)e.  Through¬ 
out  this  paper,  the  symbol  e  will  denote  a  column  vector  with  all 


its  components  equal  to  one  and  of  dimension  appropriate  to  the 
formula  in  which  it  occurs.  The  mean  service  time  y*(J)  of  the 
J-th  server  is  given  by  y ’  (J )  *  -  £(J)  S-1(J)  e. 

It  will  be  assumed  that  all  service  times  are  independent  of 
the  arrival  process.  For  any  k  *  2  successive  customer?,  the 
service  times  are  assumed  to  be  conditionally  independent,  given 
the  labels  of  the  servers  by  yhich  they  are  processed.  Under  the 
stated  assumptions,  the  model  has  a  Markov  chain >  embedded  at  the 
successive  epochs  of  arrival.  Its  states  are  of  three  types. 

For  lac,  the  state  described  by  (i,  1^,  h2,  ...,  hc)  signifies 
that,  immediately  prior  to  the  arrival,  there  are  i  customers  in  . 
the  system  and  that  for  1  s  J  sc,  the  J-th  server  is  in  the 
phase  hj  of  his  PH-distribution.  The  set  of  all  such  states  with 
a  fixed  index  i  will.be  denoted  by  I.  The  states  (c-lj  h1,  ...,  h 
are  similarly  defined,  but  one  of  the  phase  states  now  corresponds 
to  the  initial  service  phase,  selected  by  the  arriving  customer. 

The  remaining  states  correspond  to  the  case  where,  prior  to  the 
arrival,  there  are  fewer  than  c  -  1  customers  in  the  system.  The 
precise  labeling  of  these  states  is  immaterial  to  our  discussion. 

The  set  of  all  such  states  will  be  denoted  by  E.  The  states 
(i,  h1,  ...»  hc),  for  i  i  c  -  1,  1  s  J  s  c,  1  s  hj  s  v(J),  are 
listed  in  lexicographic  order.  The  states  in  E  u (c-1)  are  called 
the  boundary  states. 


In  order  to  define  the  transition  probability  matrix  P  of 
the  embedded  Markov  chain,  we  introduce  the  matrices  Pj(r,t), 
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for  r  a  0,  t  *  0,  and  1  s  J  sc.  These  matrices  satisfy  the 
differential  equations 


Cl) 


pjco.t)  -  pj co,t)  sen, 

PjCr,t)  -  PjCr,t)  S(j)  +  PjCr-l,t)  s°cn  £0), 


for  r  a  1,  with  the  initial  conditions  Pj(r,0)  *•  firQ  I,  for 
r  i  0.  Their  significance  is  the  same  as  in  C5J  or  [6]. 


We  note  that 

m 

C2)  l  P.Cr,t)  zr  -  expUSCj}  ♦  z  S°(J)  £(J)]t), 

r-0  J  “ 

for  taO,  Oszsl. 


The  matrices  Ak,  k  *  0,  are  defined  by 


(3) 


A,.  -  I 


P^r^t) 


P2 Cr2,t) 


*  Pc(rc,t)  d  F(t), 


for  k  i  0.  The  summation  is  over  all  c-tuples  (r, ,  r  ), 

t  c 

which  satisfy  *  0,  . . . ,  rQ  a  0,  r1  +  . . .  +  rc  -  k.  The  symbol 

•  stands  for  the  Kronecker  product  of  matrices.  As  shown  in  [6], 
the  matrices  Pj(r,t)  are  positive  for  r  z  1,  t  >  0,  and  1  &  J  s  c 
It  is  therefore  clear  that  the  matrices  Ak,  k  a  c,  are  positive. 

The  matrices  A^,  k  *  0,  are  of  order  m,  given  by 

OO  m  ■  n  v(J ) . 

J-l 


Prom  Formulas  (2)  and  (3),  it  follows  that 
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(5)  A*(z)  ■  Z  A.  z 

k-0  * 


|  exp{CS(l)+zS° (1)6(1) . f  exp{CS(c)+zS°(c)e(c)]t>dF(t) 


The  matrix  A*(z)  is  positive  for  0  <  z  s  1. 

The  transition  probability  matrix  P  is  now  given  by 


(6)  P 


E 

C<*1  1 

c+1 

c+2 

m 

... 

. -1 

E 

1  0 

0 

0 

0 

1 

...  j 

c-1 

!  A0 

0 

0 

0 

i 

1 

•  •  • 

w  «  «  «  « 

—  —  —  — 

—  —  —  — 

•*  —  ~  • 

~  ~  ~  i 

C 

!  ai 

Ao 

0 

0 

•  •  •  ! 

i 

c+1 

m 

1  A2 

A1 

Ao 

0 

i 

•  •  • 

c+2 

1  A^ 

A2 

A1 

Ao 

•  •  • 

sil 

!  A*» 

A3 

A2 

A1 

•  •  • 

The  elements  in  the  columns,  labeled  E  and  c-1,  are  immaterial 
to  our  discussion.  They  are,  in  general,  exceedingly  complicated 
and  depend  on  the  rule  by  which  arriving  customers  are  assigned  to 
free  servers.  Explicit,  but  highly  involved  expressions  are  given 
for  the  case  c  »  2,  in  Chapter  4  of  C8].  In  all  properly  defined 
cases,  the  matrix  P  is  irreducible. 


The 


matrix  A  »  ?  A^.  ■  A*(l),  is  a  strictly  positive,  stochastic 

k«0  K 


.ji  « •_.  j 1 a.  .. 
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matrix.  Let  0.(J )  be  the  positive  probability  vector  satisfying 

(7)  e(j)  cs(j)  +  s°(j)  i(j)]  -  o  e(j)  e  «  l, 

for  1  s  J  sc,  then  it  readily  follows  from  (5)  that  the  m-vector  w, } 
which  satisfies  v  A  ■  ir,  jr  e  »  1,  is  given  by 

(8)  jr  -  e(i)  a  ...  a  eCc). 

By  using  elementary  formulas,  proved  in  C53,  we  may  also  express 

the  vector  -  f  k  kv  e,  explicitly  in  terms  of  the  data  of  the 

k-1  K  ” 

model.  We  then  easily  verify  that 

(9)  *  B«  -  V  l  y,-1U ) . 

J-l 

As  shown  in  Chapter  1  of  [8],  the  Markov  chain  P  is  positive  re¬ 
current  if  and  only  if  1  £*  >  1>  or  equivalently 

(10)  A'"1  <  Z  v’”1U). 

J-l 

The  arrival  rate  A'”1  to  the  queue  must  be  less  than  the  combined 
service  rate  of  the  c  servers.  This  intuitive  equilibrium  con¬ 
dition  may  also  be  proved  by  applying  the  main  theorem  in  Lavenberg 
C3J. 

The  invariant  probability  vector  x  of  P  is  now  partitioned 
into  vectors  Xg,  Jq+i*  •  ••»  where  the  vectors  x^,  ii  c  -  1, 

are  m-vectors  and  the  vector  Xg  is  of  dimension  card (E) . 

It  then  follows  from  general  results,  proved  in  C 7 D  or  [8], 


that 
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(11)  x1  -  x^  r1“c+1j  for  1  i  c  -  1, 

where  the  positive  matrix  R  Is  the  minimal  nonnegative  solution 
to  the  nonlinear  matrix  equation 

(12)  R  -  ?  Rk  Av. 

k-0  k 

The  vectors  Xg  and  x^^  are  determined,  up  to  a  multiplicative 
constant,  by  solving  a  homogeneous  system  of  linear  equations. 

That  constant  is  determined  by  use  of  the  normalizing  equation 

(13)  Xg  e  +  (I-R)"1  e  »  1. 

•••••*  ..  , 

The  spectral  radius  n  *  sp(R),  is  the  unique  solution  in 
(0,1)  of  the  equation 

(14)  z  »  X(z) , 

where  X(z)  is  the  spectral  radius  of  A*(z).  The  matrix  R  is 
also  the  unique  nonnegative  solution  of  spectral  radius  less  than 
one  to  the  equation  (12). 

2.  Preliminary  Results 

The  probability  distribution  of  phase  type  with  irreducible 
representation  C£(J),  S(j)3  has  a  rational  Laplace-Stieltjes  trans¬ 
form,  given  by 

(15)  ♦j(s)  -  £(j)  CsI-SU)]"1  S°(J),  for  Re  s  4  0. 

Let  the  abscissa  of  convergence  of  $j(s)  be  -t^  <  0.  The 
function  <frj(s)  Is  then  defined,  positive  and  convex  decreasing 
on  the  interval  (-tj,«»). 
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Lemma  1. 

The  equation 

(16)  z  ♦j(s)  **  1» 

has  a  unique  real  solution  Sj  ■  tj/j(z),  for  every  z  In  (0,13. 

The  function  ’I'jC*)  satisfies  -t  <  \|ij(z)  *  0»  and  Is  strictly 
Increasing  on  (0,13.  Moreover  ^(0+)  ■  -tj ,  and  tpj  (1)  «  v*'"1(J). 

The  quantity  ^(.z)  is  the  eigenvalue  of  maximal  real  part 
of  the  matrix  S(J)  +  zS° (J )  jKj).  The  corresponding  left  eigen¬ 
vector  u(J,z),  normalized  by  u(j,z)  e  *  1,  Is  given  by 

(17)  u ( J , z )  -  z(z-l)“1  ^(z)  0(j)C«J(z)I-S(j)3“1,  for  0  <  z  <  1, 

*  0(J),  for  z  *  l. 


Proof 

Essentially  the  same  results  were  proved  in  C113.  For  ease 
of  reference,  we  repeat  the  proof.  The  first  set  of  properties  of 
#j(z)  follow  readily  from  consideration  of  the  graph  of  $j(s) 
and  from  the  equation  (16).  The  equation 

(18)  u(j,z)  CS(J)  +  z  S°(J)  £(J)]  -  V»j(z)  u(J,z), 

leads  to 

u(J,z)  ■  z  Cu(J,z)  S 0 ( J ) 3  £(J)  Cipj  (z)I  -  S(J)3“1. 

Postmultiplication  by  S°(J)  leads  to  z4>j[t|/j(z)3-l.  The 
inner  product  u(J,z)  S°(J)  does  not  vanish,  since  the  matrix 
i(ij(z)  I  -  S ( J )  is  nonsingular  for  0  <  z  t  1. 
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The  vector  u(J,z)  =  J  expC-4»j (z)t ]*£(J )  exp[S(J)t3  dt,  is 
positive,  since  the  vector0  £(J )  exp[S(j)t]  is  positive  for  t  >  0, 
as  was  shown  in  [63.  This  implies  that  the  eigenvalue  4».(z)  of 

J 

the  irreducible  stable  matrix  S(J)  +  z  S°(j)  g(j),  is  the  eigen¬ 
value  of  maximal  real  part.  The  normalization  u(J,z)  e  *  1, 
readily  yields  (17). 

Lemma  2 

The  maximal  eigenvalue  X(z)  of  A*(z)  is  given  by 

c 

(19)  X(z)  «  f[  _  Z  p,(z)],  for  0  <  z  s  1, 

J-l  J 

where  f(*)  is  the  Laplace-Stieltjes  transform  of  the  interarrival 
time  distribution  F(-).  The  corresponding  eigenvector  u(z)  is 
given  by 

(20)  u(z)  ■  u(l,z)  a  ...  a  u(c,z). 

Proof 

The  vector  u(z)  is  clearly  positive  and  satisfies  u(z)  e  ■  1, 
It  readily  follows  from  (18)  that 

u(J,z)  exp  {[S(j)  +  z  S°(J)  jg.(  J )  3t  >  ■  exp  C^j(z)t3  u(J,z), 

and  hence  by  (5),  that 

c 

u(z)  A*(z)  ■  f  C  -  I  i|>.(z)3  u(z). 

J-l  3 

This  clearly  implies  (19)  and  completes  the  proof. 

Let  now  n  be  the  unique  solution  in  (0,1)  of  the  equation 
(14).  The  vector  u(ti)  ■  u(l,n)  •  ...  ®  u(c,n),  Is  then  given  by 
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u(j,n)  «  ntn-i)”1  ^(n)  i(j)C^Cn)  I  -  SQ)]"1,  for  1  s  j  sc. 

As  shown  In  Chapter  1  of  [8],  the  vector  u(n)  is  also  the  left 
eigenvector  of  the  matrix  R,  corresponding  to  its  Perron  eigen¬ 
value  n . 


3.  Asymptotic  Behavior  of  the  Queue  Length  Density 
Theorem  1 

The  stationary  density  of  the  queue  length  at  arrivals  satisfies 

(21)  I  x-  e  *  (l-rO^Cx,  ,z)  nk"c+1  +  o(nk),  as  k  ♦  », 

i«k  -1  “  -c~i 

where  z  is  the  right  eigenvector  of  R,  corresponding  to  the 
eigenvalue  n  and  satisfies  u(n)  z  *  1. 


Proof 

Let  us  write  u  for  u(r\).  A  classical  property  of  irre¬ 
ducible,  nonnegative  matrices  now  yields  that 

R1  *  n1  z.  u  +  as  i 

Since 

40 

I  *4  e  -  x  ,  Rk"c+1(I-R)_1  e, 
i-k  ~ 

Formula  (21)  readily  follows. 


Remark 

We  clearly  also  have 

**  *  Sq.i  R1"0*1  *  n1“c+1  £ )  y.  +  otn1),  as  l  +  ®, 


(22) 


This  implies  that 
(23)  5t 

.  u ,  as  i 

-i  *  " 

The  left  hand  side  is  the  conditional  probability  density  of 
the  c-tuple  of  service  phases,  given  that  the  arriving  customer 
finds  i  customers  in  the  queue.  For  large  i,  we  find  that  the 
Joint  conditional  contribution  of  the  c  residual  service  times 
has  the  limit  distribution  with  c-fold  Laplace-Stieltjes  transform 

h(s1,...,sc)  -  }HCs1I-S(l)3“1S0Cl)8...f[scI-S(c)]"1S0(c)> 

-  n  u(j,n)  Cs.I-SU)]"1  S°(J). 

J-l  ”  J 

We  see  that  for  large  i,  the  conditional  residual  service  times 
in  the  c  servers,  given  that  the  arriving  customer  finds  i 
customers  in  the  system,  are  approximately  Independent .  The 
marginal  distribution  of  the  residual  service  time  for  the  J-th 
server  is  then  approximately  the  PH-distribution  with  representa¬ 
tion  Cu(J,n),  S(J)],  for  1  s  J  sc. 

4.  Asymptotic  Behavior  of  the  Stationary  Waiting  Time  Distribution 
at  Arrivals. 

Once  the  vector  x  .  and  the  matrix  R  are  known,  the 
probability  distribution  W(*)  of  the  waiting  time  at  arrivals 
under  the  first-come,  first-served  discipline,  may  be' computed  by 
solving  a  finite,  highly  structured  system  of  differential  equations 
with  constant  eoefflcients.  This  will  be  shown  by  generalizing  the 
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proof,  given  for  the  case  of  the  GI/PH/1  queue,  in  C8)  or  C93» 
For  multi-server  queues,  the  practical  utility  of  this  result  is 
severely  limited  by  the  dimension  m.  The  proof  of  this  result  is 
however  essential  to  the  derivation  of  the  desired  asymptotic 
formula  for  W(*)* 

Let  the  matrices  C  and  D  be  defined  by 
C  »  S(l)  9  S(2)  ®  ...  9  S(c ) , 

(24) 

D  -  S°(l)  £(1)  9  S°(2)  6(2)  «  ...  «  S°(c)  6(c), 


where  the  sumbol  9  is  the  Kronecker  sum  of  matrices  [ID.  The 
distribution  W( • )  may  be  viewed  as  the  distribution  of  the  time 
till  absorption  into  the  state  *  in  the  Markov  process  with 
generator  Q,  given  by 


(25)  Q 


in  which  the  initial  probability  vector  is  given  by 


2E  £  +  Jci  S.,  x^  R,  R‘ 


AHUM  IWW 
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It  is  readily  seen  that  the  Laplace-Stieltjes  transform  w(s) 
of  W(*)  is  given  by 

(26)  w(s)  -  x-  e  +  x  E  R1  C  (sI-O^D]1  e 

i«0  ~ 

-  2e  e  +  x^  ?»(s)  e, 

where  f*(s)  is  the  square  matrix  of  order  m,  which  satisfies 

(27)  ¥»(a)  -  I  +  R  ¥»(s)  (sI-C)”1  D. 

The  equation  (27)  is  now  transformed  in  the  same  manner  as  discussed 
in  191 •  If  ¥(«)  Is  the  matrix  of  mass  functions  with  Laplace- 

p 

Stieltjes  transform  ¥»(s)  and  £(x)  is  the  m  -vector  obtained  by 
forming  the  direct  sum  of  rows  of  ¥(x),  then  we  derive  from  (27) 
that 

(28)  ♦(*)  -  v  -  v(I|C  +  RTf  D)_1(I*I  -  expC (X9C  +  RT»D)x]>(RT«D) , 

2 

for  x  2  0.  The  vector  v  is  the  m  -vector  obtained  by  forming 

T 

the  direct  sum  of  the  identity  matrix.  R  is  the  transpose  of  the 
matrix  R. 

We  now  set  v°  -  -  v(I*c  +  R1*!)"1,  and  e(x)  - 
v°  expC(ItC  +  R^tD)x3,  for  x  a  0.  The  m  x  m  matrices  V°  and 
0(x)  have  the  vectors  v°  and  0(x)  as  the  direct  suras  of  their 
respective  rows.  They  satisfy  the  equations 

(29)  V°  C  ♦  R  V°  D  ■  -  I, 
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and 

(30)  e*(x)  -  e(x)  c  +  r  e(x>  d,  e(o>  -  v°, 

for  x  *  0 . 

By  virtue  of  (28),  the  matrix  ¥(x)  is  then  given  by 

(31)  ¥(x)  -  I  ♦  RV°  D  -  R  0(x)  D,  for  x  *  0. 

The  distribution  W(*)  is  given  by 

(32)  W(x)  -  Xg  e  +  x e-1  e  +  xc-1  R  V°  D  e  -  x^  R  0(x)  D  e, 

for  x  a  0.  This  expression  may  be  further  simplified.  We  post- 
multiply  in  (29)  by  e  and  note  that  C  e  +  D  e  *  0 .  This  yields 
that  V°  D  e  •  (I-R)"1  e.  Upon  substitution  into  (32),  we  readily 
obtain 

(33)  W(x)  -  1  -  x^  R  ©(x)  D  e,  for  x  *  0. 

We  see  that  the  probability  distribution  W(*)  may,  in  principle, 
be  computed  by  first  evaluating  the  matrix  V°  and  then  solving 
the  matrix-differential  equation  (30).  In  order  to  obtain  the 
asymptotic  formula,  we  neod  a  number  of  preliminary  lemmas.- 

Lemma  3 

The  matrix  C  +  n  D,  given  by 

(34)  CMD*  CS(1)  +  n  S°(l)  £(1)3«  ...  «CS(c)  +  n  S°(c)  £(c)], 

is  an  Irreducible,  stable  matrix.  Its  eigenvalue  -  ?  of  maximal 
real  part  is  given  by 

c 

£  •  I  U>,(n). 

J-l  J 


(35) 
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The  corresponding  left  eigenvector  Is  given  by  u  ■  u(n).  The 
corresponding  right  eigenvector  u°,  normalized  by  u  u°  ■  1,  Is 
given  by  the  Kronecker  product  u°  ■  u°(l)  ft  ...  t  u°(c),  where 


C36)  u°U) 


_  n-1 . . 

(n)  *£(j  )c^  (n)i-s(j  )]”2s°(j ) 


C*j  (nJl-sOJl’Vu), 


for  1  s  J  sc. 


Proof 

Since  each  of  the  matrices  SCJ)  +  n  S°(J)  £(J),  1  s  J  sc.  Is 
an  Irreducible  stable  matrix,  so  Is  the  matrix  C  +  r\  D,  [111.  The 
matrix  C  +  n  D  Is  the  sum  of  c  matrices  of  the  form 


If...  gCS(J)  +  n  S°(J)  *  jgCj )]  I  I  I  ...  •  I. 

This  readily  yields,  by  (18),  that 

c 

u(C  +  n  D)  •  L  *  u  . 

J-l  J 

The  vector  u°(J)  Is  clearly  a  right  eigenvector  of 
S(J)  +  n  S°(J)  *  Jl<j),  corresponding  to  Furthermore 

o  a  0 

u(4,n)  u  CJ)  *  1*  Since  u  u°  »  n  u(J,n)  u°(J)  ■  1,  the  proof 

~  J-l  “ 

Is  complete. 


Lemma  M 

The  matrix  I  •  C  +  •  D  has  nonnegative  off-diagonal 

elements  and  is  Irreducible.  Its  eigenvalue  of  maximal  real  part 
Is  -  £.  The  corresponding  left  and  right  eigenvectors  are  re¬ 
spectively  given  by  z,T  f  u  and  uT  •  u°.  Their  Inner  product 
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Is  one. 

Proof 

T 

The  off-diagonal  elements  of  I  •  C  +  R  •  D  are  clearly 

T 

nonnegative.  The  irreducibility  of  the  matrix  I  £  C  +  R  §  D 
follows  from  the  positivity  of  R  and  the  Irreducibility  of  the 
representations  of  the  service  time  distributions. 

We  have 

(zT  t  u)Cl  C  C  +  RT  ft  D)  -  zT  ft  u  C  +  ft  zT  ft  u  D  -  zT  C  u(C  +  n  D) 

-  C(zT  8  u), 

and  a  similar  calculation  for  the  right  eigenvector.  Since  both 
eigenvectors  are  positive,  -  C  Is  the  eigenvalue  of  maximal  real 
part  [2], 

Lemma  5 
The 

(37) 

Also 

(38) 

Proof 

It  follows  from  the  definition  of  v°  that 


vector  v°  ■  -  v(I  C  C  +  RT  ®  D)”1,  satisfies 
v°(uT « u°)  -  r1. 

u  D  e  «  5(1  -  n)*1. 


However 


v(uT  0  u°)  -  eg,  e^] 


U1  E 


u2  — 


u_  u 

Til  — 


I  %(£».  E°)  ■  E  E°  “  1< 

v«i  v 


The  vectors  e  are  the  m  unit-vectors  of  dimension  m.  This 

— V 

proves  Formula  (37 ) - 

A  typical  term  of  D  e  is  the  Kronecker  product 
e  0e  0...  0S°(J)0e  ...  0  e .  Premultiplication  by  u  -  u(l)  0  . 
0u(c)  yields 

u(j )  S°CJ)  -  n(n  -  l)"1  (n)  £(J)C^Cn)l  -  S(J)]"1  S°Q)  - 


(n  -  l)-1i|>j(n) 


so  that 


u  D  e  -  (u  -  l)"1  t  *,(n)  »  5(1  -  n)"1. 

>1  J 


Theorem  2 

The  waiting  time  distribution  W(*)  satisfies 
(39)  1  -  W(x)  -  nCl  -  h)"1  Cxc-1  z)  e“5x  +  o(e"?x),  as  x 

Proof 


It  follows  from  the  definition  of  the  vector  8(x)  and  the 

T 

properties  of  the  matrix  I  0  C  +  R  0  D,  that 
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8(x)  -  v°C(uT  >  u°)  •  (zT  ft  u)]e"5x  +  o(e“5x) 

■  •  u)e~*x  +  o(e"^x),  as  x  •». 

T 

The  vector  £  •  u  is  the  direct  sum  of  the  rows  of  the  matrix 
z  •  u.  The  preceding  formula  may  therefore  be  equivalently  written 
as 

©(x)  *  C^Cz  •  u)e"^x  +  o(e~^x),  as  x  -*■  <*. 
Substitution  into  (33)  yields 
1  -  W(x)  ■  5”1Cxc_1  R  z)(u  D  e)e“^x  +  o(e"^x),  as  x  -*•  ». 

Since  R  z  *  t)  z,  and  by  using  Formula  (38),  we  readily  obtain  (39). 

Remark 

Results,  similar  to  those  in  Theorems  1  and  2,  may  be  proved 
for  the  stationary  distributions  of  the  queue  length  and  the 
waiting  time  at  an  arbitrary  time.  The  proofs  proceed  along  the 
same  lines  as  in  the  single  server  case,  discussed  in  C93.  The 
same  decay  parameters  n  and  g  are  obtained,  but  the  multi¬ 
plicative  constants  are  different, 

5.  Computational  Procedure  and  Applications 

The  decay  parameters  n  and  g  may  be  computed  together  by 
elementary  algorithms.  There  are  various  alternative  methods.  It 


Is  advislble  to  solve  the  equation 

(40)  z  -  f£-  ^  fy(z)J, 

for  n  In  (0,1)  by  a  method  which  does  not  Involve  derivatives. 

The  secant  or  bisection  methods  may  be  Implemented  with  equal  ease. 

At  each  stage  of  the  computation,  we  have  two  values  and 

z2  satisfying 

0  <  {-  X  ♦,<*!>]  «  '  ‘  {-  X  W]  ‘  22  '  11 

since  the  right  hand  side  of  (40)  Is  Increasing.  As  the  next  trial 
value  z’  Is  obtained,  either  by  bisection  or  the  secant  method, 
the  corresponding  values  ifij(z’),  1  s  J  s  c,  are  computed  by 
solving  the  equations 

(41)  z  £(J  ) C (z ) I  -  S(J)]”1  S°(J)  -  1, 

for  their  unique  solutions  In  the  Intervals  (-Tj,0),  1  s  J  sc. 

One  clearly  only  solves  those  equations  which  are  actually 
different.  The  monotonicity  properties  of  the  ^j(z),  proved  In 
Lemma  1  are  useful  In  solving  the  equations  (4l).  When  the 
interval  (z1,  z2),  which  brackets  n  Is  sufficiently  small,  we 
evaluate  a  final  value  n,  which  Is  the  computed  value  of  r>. 

A  ®  A 

The  computed  value  C  of  £  is  obtained  by  setting  £  —  Z  ^.(n). 

J»1  J 

For  many  PH-dlstrlbutions  of  interest,  the  Laplace-Stieltjes 
transform  la,  of  course,  explicitly  available,  so  that  the 
equations  (4l)  can  then  be  written  In  a  computationally  more 
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convenient  form. 

Except  for  queues  with  a  very  small  number  of  servers  and 
then  only  for  PH-distributions  with  few  phases,  the  computation  of 
the  matrix  R,  and  hence  of  the  vectors  5<J^1  and  z,  is  not 
practically  feasible.  Even  without  explicit  knowledge  of  the 
constant  x^_^  z_t  the  asymptotic  results  of  Theorems  1  and  2  have 
practical  uses. 

With  n  and  £  so  easily  computable,  these  results  may  be 
used  to  test  the  merits  of  simulation  procedures  for  the  queue 
length  and  waiting  times  in  multiserver  queues.  The  estimates  of 
-log  Cl  -  W(x)],  for  example,  should  for  large  x  lie  approximately 
on  a  straight  line  of  slope  £.  Assuming  that  the  simulation  pro¬ 
cedure  can  correctly  identify  the  parameter  £,  it  should  also  be 
sufficiently  accurate  to  provide  a  good  estimate  for  the  intercept 
of  the  linear  asymptote  of  -log  Cl  -  W(x)3.  We  will  then  have  an 
estimate  of  x^^  z»  which  may  be  used  in  the  asymptotic  formulas 
to  provide  estimates  of  tall  probabilities  for  the  queue  length  and 
waiting  time. 

As  a  point  of  theoretical  interest,  it  appears  likely  that  the 
asymptotic  results  of  Theorems  1  and  2  remain  valid  for  the  OI/Q/c 
queue  with  heterogeneous  servers,  provided  that  each  of  the  c 
service  time  distributions  have  a  Laplace-Stieltjes  transform  with 
a  negative  abscissa  of  convergence.  This  may  probably  be  proved  by 
appropriate  continuity  arguments  and  the  approximation  of  the 
service  time  distributions  by  PK-dlstrlbutlons.  This  matter,  as 


well  as  the  applications  to  simulation  methodology,  will  be  taken 
up  elsewhere. 
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